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Results  of  HEMP  code  calculations  of  simple  problems  are  compared 
with  exact  results  to  determine  the  accuracy  with  which  this  numerical 
technique  solves  initial -boundary  value  problems.  More  complicated  prob 
lems  of  metal  liners  accelerated  by  explosives  were  also  calculated  with 
this  code,  and  these  results  were  compared  with  experimental  results  to 
determine  how  well  the  code  simulates  these  physical  processes.  It  is 
concluded  that  HEMP  calculations  are  accurate  enough  to  be  useful  as  a 
tool  for  predicting  the  gross  motion  of  metals  accelerated  by  detonating 
explosives. 
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I.  INTRODUCTION 


The  HEMP1*  code  is  a  computer  program  which  is  formulated  to  ap¬ 
proximately  solve  the  conservation  equations  of  mass,  momentum,  and  en¬ 
ergy  subject  to  appropriate  initial  and  boundary  conditions.  The  ma¬ 
terials  under  consideration  are  described  as  hydrodynamic  materials 
(fluids)  or  elastic-perfectly  plastic  materials.  An  energy  release 
routine  which  simulates  the  detonation  of  a  high  explosive  is-  also  in¬ 
cluded.  This  code  contains  the  essential  features  necessary  to  solve 
explosive-metal  acceleration  problems  which  are  relevant  to  the  design 
of  warheads.  Since  we  plan  to  use  this  code  to  study  and  predict  the 
performance  of  warheads,  we  first  want  to  gain  an  understanding  of  the 
accuracy  obtainable  from  it.  Two  important  questions  regarding  the 
accuracy  of  finite-difference  codes  are:  (1)  how  well  does  the  numer¬ 
ical  scheme  solve  the  governing  system  of  partial  differential  equa¬ 
tions,  and  (2)  how  well  do  these  equations  model  the  physical  process? 
To  obtain  some  sort  of  answer  to  these  questions,  we  first  compared 
numerical  solutions  with  exact  solutions  to  simple  problems,  and  then, 
for  more  complicated  problems,  we  compared  numerical  solutions  with 
experimental  results. 

In  Section  II,  a  sequence  of  simple  problems  which  possess  exact 
solutions  is  solved  numerically  with  the  code,  and  comparisons  between 
the  exact  and  approximate  solutions  are  made.  Problems  Involving  per¬ 
fect  fluids,  elastic  solids,  elastic-plastic  solids,  and  detonating 
explosives  were  considered. 

In  Section  III,  comparisons  between  experimental  test  results  and 
HEMP  code  calculations  are  made  for  explosive-metal  systems.  Conclu¬ 
sions  from  these  comparisons  are  presented  in  Section  IV.  Although  the 

*Referenoea  are  listed  on  page  35. 
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problem  formulation  and  method  of  solution  used  in  the  HEMP  code  have 
been  documented  by  the  originators1*2  ,  a  very  brief  outline  of  the 
formulation  is  presented,  for  completeness,  in  Appendix  A. 


II.  HEMP  SOLUTIONS  COMPARED  WITH  EXACT  RESULTS 

In  this  section,  the  results  of  HEMP  calculations  applied  to  a 
series  of  simple  problems  are  presented.  For  these  problems,  the  exact 
solutions  are  well  known  and  can  easily  be  evaluated.  These  exact  solu¬ 
tions  are  then  compared  with  the  HEMP  results  in  order  to  determine  how 
well  the  governing  equations  are  satisfied  by  the  HEMP  solutions. 

As  a  first  illustration,  the  problem  of  a  high  pressure,  ideal 
gas  initially  trapped  in  a  tube  and  then  released  was  calculated  with 
three  different  mesh  sizes.  The  results  are  illustrated  in  Figure  1. 
For  this  and  subsequent  figures,  the  notation  (P)  pressure,  (p)  density, 
(e)  specific  internal  energy,  (y)  a  constant,  (G)  shear  modulus,  (ay) 
yield  strength  in  tension,  and  (Vo)  initial  velocity  is  used.  The  sol¬ 
id  line  on  the  pressure-position  graph  indicates  the  exact  solution  at 
10  ysec  after  release  of  the  gas,  and  the  three  numerical  results  are 
indicated  by  the  points.  As  the  number  of  points  in  the  numerical  mesh 
is  increased  the  numerical  solution  approaches  the  exact  solution,  and 
convergence  of  the  numerical  scheme  is  indicated.  All  solutions  repre¬ 
sent  the  plane  rarefaction  wave  reasonably  well.  The  variation  of  prop¬ 
erties  across  the  tube,  for  the  HEMP  results,  was  negligible.  For  com¬ 
putational  purposes,  this  one-dimensional  problem  can  be  generalized 
slightly  by  using  a  mesh  which  Is  initially  skewed  as  shown  in  Figure  2. 
The  plane  rarefaction  then  passes  obliquely  across  the  zones,  and,  for 
computational  purposes,  the  problem  is  two-dimensional.  The  numerical 
solution  is  degraded  only  slightly  because  of  this  two-dimensional  ef¬ 
fect,  as  Indicated  by  the  comparison  between  calculations  with  the  rec¬ 
tangular  and  non-rectangular  meshes  shown  in  Figure  2.  The  most 
noticable  result  of  the  non-rectangular  zoning  appears  in  the  position 
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Figure  1.  Numerical  solution  for  the  one-d1mt>is1onal  expansion  of  a 
high  pressure  gas  using  three  different  finite-difference 
zonlngs 
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of  the  gas  boundary.  However,  for  practical  purposes,  HEMP  solutions 
appear  to  be  reasonably  accurate  for  expansion  waves  in  an  ideal  gas. 

Since  an  interest  in  the  flow  of  metals  also  exists,  a  series  of 
one-dimensional  (uniaxial  strain)  impact  problems  was  solved  with  the 
HEMP  code.  Figure  3  shows  a  typical  Hugoniot  curve  (solid  line)  for 
an  elastic-perfectly  plastic  material3.  The  dotted  line  indicates  the 
straight  line  extension  of  the  elastic  portion  of  the  Hugoniot.  Three 
different  wave  structures,  depending  on  the  magnitude  of  the  stress 
jump,  can  exist.  If  a  stress  jump  of  ci  is  present,  where  ctj  is  below <rv 
the  yield  point  stress,  an  elastic  wave  with  a  single  stress  jump  will 
be  propagated.  A  double  wave  structure,  elastic  precursor  followed  by 
a  plastic  shock,  will  propagate  if  the  stress  jump  is  a2,  where  a2  is 
a  stress  above  the  yield  point  value  but  below  the  intersection  of  the 
dotted  line  with  the  Hugoniot.  A  stress  jump  a3  above  the  intersection 
point  will  again  result  in  a  single  stress  jump  (shock). 

The  numerically  calculated  pressure  distribution  resulting  from 
the  elastic  Impact  of  two  identical  blocks  Is  shown  in  Figure  4  along 
with  the  exact  solution.  The  zoning  and  initial  conditions  of  the 
problem  are  also  shown.  Boundary  conditions  of  no  vertical  velocity  on 
the  upper  or  lower  surfaces  were  used  in  the  numerical  solution.  The 
HEMP  values  oscillate  somewhat  (±6%)  about  the  exact  value  of  8  Kb. 

A  smoother  solution  to  this  impact  problem  may  be  obtained  by  Increas¬ 
ing  the  number  of  zones  and  by  inserting  a  transition  zone  to  smooth 
the  Initial  velocity  discontinuity.  Results  of  calculations  with 
these  changes  are  shown  In  Figure  5  along  with  the  initial  velocity 
distribution.  The  improvement  in  the  quality  of  the  solution,  compared 

to  the  solution  shown  in  Figure  4,  is  apparent.  This  problem  may  also  be 

treated  as  a  two-dimensional  problem  by  introducing  the  non-rectangular 

zoning  indicated  in  Figure  6.  The  degradation  of  the  solution  due  to  two- 

dimensional  effects  can  be  seen  by  comparing  Figure  6  with  Figure  4.  The 

oscillation  In  the  numerical  solution  for  the  non-rectangular  zoning  case 

has  increased  in  amplitude  by  a  factor  of  two  over  the  one-dimensional  case 

However,  the  solution  can  be  made  smoother  by  using  more  zones.  A 
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Figcre  3.  Schematic  of  a  Hugoniot  Curve  for 
Elastic  -  Plastic  Material 
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Figure  4.  The  Initial  configuration  and  pressure  distribution  3  psec 
after  the  one-dimensional  Impact  of  elastic  blocks 
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Figure  5.  Improvement  in  the  numerical  solution  of  the  elastic 
impact  problem  by  using  an  initially  fine  zoning  and 
a  transition  zone 
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Figure  6.  Effect  of  Initial  zone  shape  upon  numerical  solution  of 

the  elastic  Impact  problem 
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numerical  solution  using  40  axial  zones  in  the  target  material  instead 
of  20  is  shown  in  Figure  7  where  both  the  pressure  and  particle  velo¬ 
city  distribution  are  indicated.  It  has  been  suggested  that  oscilla¬ 
tions  in  these  solutions  could  be  additionally  reduced  by  introducing 
a  linear  artificial  viscosity  term1. 

If  tie  impact  velocity  is  sufficiently  increased,  the  two  wave 
structure,  indicated  by  the  stress  level  of  Figure  3,  can  be  formed 
in  the  Impacting  materials.  Figure  8  illustrates  a  part  of  the  numer¬ 
ical  solution  for  this  case.  If  the  impact  velocity  is  further  in¬ 
creased,  the  single  wave  structure,  corresponding  to  the  stress  level 
<r3,  may  be  formed.  The  numerical  results  for  this  high  velocity  Im¬ 
pact  case  are  shown  in  Figure  9. 

The  conclusion  from  this  series  of  Impact  problems  Is  that  the 
H£MP  results  do  represent,  within  reasonable  accuracy,  the  steady-state 
wave  profiles  which  occur  with  the  elastlc-perfectly  plastic  material. 

To  simulate  the  detonation  of  a  high  explosive,  an  "explosive 
burn  routine"  is  Included  In  the  code.  The  chemical  energy  released  in 
the  detonation  process  Is  stored  In  each  zone  as  an  Initial  Internal 
energy.  Release  of  this  energy  Is  started,  In  each  zone,  when  the  cal¬ 
culated  position  of  the  detonation  reaches  the  zone  center.  Tu  test 
this  portion  of  the  code,  the  problem  of  a  one-dimensional  Chapman- 
Jouquet  (C-J)  detonation  wave  propagating  down  a  closed-end  tube  was 

numerically  solved.  Figure  10  shows  the  exact  pressure  distribution 
compared  with  two  numerical  solutions  (100  zones  and  20  zones).  Both 
solutions  fall  to  produce  the  peak  pressure  (C-J  pressure).  The  peak 
pressure  of  the  finely  zoned  solution  Is  about  10X  low.  For  this  prob¬ 
lem,  the  left  boundary  is  fixed.  If  the  left  boundary  Is  moved  as  a 
piston  to  the  right  with  the  C-J  particle  velocity,  the  peak  pressure 
of  the  numerical  solution  does  reproduce  the  exact  C-J  value  for  both 
zonings. 
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Figure  7.  Improvement  In  the  numerical  solution  of  the  elastic  problem 
with  skewed  zones  by  refining  the  zoning 
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III.  HEMP  SOLUTIONS  COMPARED  WITH  EXPERIMENTAL  RESULTS 

The  previous  section  demonstrated,  for  a  few  cases,  the  accuracy 
with  which  the  code  solves  the  partial  differential  equations  and  re¬ 
lated  jump  conditions.  In  this  section,  the  code  calculations  are  com¬ 
pared  with  experimental  results  to  determine  the  overall  accuracy  with 
which  the  code  simulates  the  physical  process  of  metal  acceleration  by 
explosive  loading. 

The  problem  of  a  cylindrical  charge  of  explosive  (Composition  B) 
accelerating  an  aluminum  disc  was  chosen  for  the  first  comparison.  The 
experimental  configuration  as  well  is  the  experimental  and  numerical 
results  are  shown  In  Figure  11,  Here,  the  motion  of  the  free  surface 
of  the  aluminum  disc  was  measured  as  a  function  of  time  with  a  streak 
camera.  The  calculations  were  based  on  the  JWL  equation  of  state4  for 
the  explosive  and  an  equation  of  state  described  by  Wilkins1  for  the 
aluminum  disc.  The  calculated  results  agree  reasonably  well  with  the 
experimental  values  (6  1/2%  maximum  error). 

Figure  12  Illustrates  the  results  of  a  test  which  was  designed 
to  measure  the  terminal  velocity  of  preformed  steel  fragments  projected 
from  the  end  of  a  Composition  B  charge.  The  experimental  velocities 
were  reported  by  Taylor5.  The  calculations  were  based  on  the  JWL  equa¬ 
tion  of  state  for  Composition  B  and  the  equation  of  state  which  Is  in¬ 
corporated  In  the  code  for  Iron.  To  model  the  preformed  fragments,  a 
hydrodynamic  description  (yield  strength  equal  to  zero)  was  used  with 
the  minimum  pressure  set  at  zero.  The  HEMP  results  duplicate  the  ex¬ 
perimental  values  within  about  6%.  If  one  would  use  a  Gurney  type 
analysis  to  predict  velocltes  for  this  configuration,6  the  predicted 
velocity  would  be  about  2.8  mm/ysec,  or  a  minimum  error  of  about  32%. 

A  third  comparison  between  calculated  and  experimental  results  Is 
Illustrated  In  Figures  13  and  14.  A  steel  cylinder  filled  with  Com- 
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Figure  11.  Calculated  and  experimental  results  for  the  accelera¬ 
tion  of  an  aluminum  disc  propelled  by  a  composition  B 
explosive  charge 
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Calculated  and  experimental  results  for  the  accelera¬ 
tion  of  preformed  fragments  from  the  end  of  an  explo¬ 
sive  charge 
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Figure  13.  Calculated  and  experimental  results  for  the  speed 
distribution  of  fragments  propelled  from  an  explo 
sively  loaded  cylinder 
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Figure  14,  Calculated  and  experimental  results  for  projection 
angle  distribution  of  fragments  propelled  from  an 
explosively  loaded  cylinder 
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position  8  is  detonated  at  one  end  and  the  resulting  fragment  velocities 
are  both  measured  and  calculated.  In  the  calculations,  the  steel  is 
considered  as  an  elastic-perfectly  plastic  solid  up  to  the  point  of  ex¬ 
pansion  where  the  pressure  becomes  less  than  a  predetermined  mlnlmun 
value.  From  that  point  on,  the  steel  Is  treated  as  a  hydrodynamic  ma¬ 
terial  in  an  attempt  to  model  the  fragmented  liner.  The  calculations 
were  stopped  when  It  was  observed  that  the  velocities  of  the  steel  lin¬ 
er  had  stabilized  (about  45  psec  after  Initiation}.  The  experimental 
fragment  velocities  were  measured  from  fragment  positions  obtained  at 
300  psec  and  at  600  usee.  These  velocities  are  compared  In  Figure  13 
as  a  function  of  the  Initial  position  along  the  cylinder.  Figure  14 
compares  the  calculated  and  observed  fragment  projection  angles.  Fig¬ 
ures  13  and  14  indicate  good  agreement  between  calculated  and  ex¬ 
perimental  results. 

Figures  15,  16  and  17  show  the  final  comparison  between  cal¬ 
culated  and  experimental  results.  The  top  portions  of  these  photographs 
show  flash  radiographs  of  the  event  while  the  lower  portions  show  the 
calculated  configurations.  The  explosive-metal  system  consists  of  a 
spherical  cap  of  mild  steel  in  contact  with  a  cylindrical  explosive  (Com¬ 
position  B)  and  a  booster.  The  metal  is  accelerated  from  the  left  face 
of  the  explosive  which  is  initiated  on  the  right  end  of  the  booster. 
Figures  15,  16  and  17  show  various  stages  of  the  acceleration  pro¬ 
cess  to  the  point  where  the  numerical  solution  could  not  be  continued 
without  rezoning.  In  these  figures,  the  relative  horizontal  location  of 
both  the  radiographs  and  the  calculated  grid  is  correct.  The  last 
comparison  (31.8psec)  shows  the  configuration  after  the  liner  has  moved 
about  one  and  a  half  calibers.  The  difference  between  the  calculated 
velocity  at  this  time  and  velocities  estimated  from  the  radiographs  is 
6  72*.  For  this  series  of  calculations,  the  JWL  equation  of  state  was 
used  to  describe  the  explosive  products.4  The  liner  material  was  model¬ 
ed  as  an  elastic-perfectly  plastic  material  with  an  equation  of  a  state 
of  the  Mie-Griinelsen  form.  To  simulate  the  edge  spall  which  occurs  in 
the  lifter,  a  failure  criterion  based  upon  pressure  was  used.  If  the 
pressure  in  a  zone  goes  below  a  minimum  value  (-10Kb),  the  material  in 
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Figure  15.  Calculated  and  experimental  configuration  of  a  steel 
spherical  cap  accelerated  by  an  incontact  explosive 
at  the  Instant  of  initiation  (right)  and  7  usee  after 
Initiation  (left) 
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Figure  17.  Calculated  and  experimental  configuration  of  a  steel 
spherical  cap  accelerated  by  an  incontact  explosive 
at  22.5  ysec  (right)  and  31.8  ysec  (left)  after  in¬ 
itiation 
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that  zone  Is  considered  failed,  and  the  yield  strength  and  shear  modulus 
are  set  to  zero  In  that  zone  for  the  remainder  of  the  calculation. 

Also,  the  subsequent  pressure  In  that  zone  Is  not  allowed  to  become  neg¬ 
ative.  At  12. 5  usee,  spalling  of  the  edges  of  the  liner  Is  apparent  In 
the  radiograph.  It  Is  also  Indicated  In  the  calculation  by  the  severely 
distorted  zones.  In  the  calculations,  all  materials  are  forced  to  re¬ 
main  continuous;  therefore,  spalled  regions  are  Indicated  by  grossly 
distorted  zones,  see  Figure  17.  In  order  to  continue  these  calcula¬ 
tions  to  31  psec.  It  was  necessary  to  disregard  the  time  step  calculation 
for  the  spalled  zones.  The  calculated  motion  of  the  remainder  of  the 
grid  should  be  valid. 


IV.  SUMMARY  AND  CONCLUSIONS 

In  an  attempt  to  gain  an  understanding  of  the  accuracy  obtained 
from  the  HEMP  code,  five  problems  were  solved  numerically  with  It,  and 
these  solutions  were  compared  with  the  exact  solutions  to  the  same  prob¬ 
lems.  These  problems  were  selected  because  they  exercised  the  subrou¬ 
tines  which  are  necessary  to  solve  problems  of  metal  liners  accelerated 
by  high  explosives.  Based  upon  this  comparison  the  numerical  technique 
appears  to  generate  solutions  which  agree  reasonably  with  the  exact 
solutions.  Problems  which  are  Inherently  two-dimensional  require  a 
finer  zoning  than  one-dimensional  problems.  If  the  same  quality  of 
solution  Is  desired.  Generally,  If  a  coarse  zoning  Is  used,  structure 
which  appears  In  the  exact  solution  tends  to  be  smoothed  and  numerical 
values  may  tend  to  oscillate  about  the  correct  values.  Also,  solutions 
to  problems  which  are  essentially  hydrodynamic  flows  appear  to  be  better 
than  solutions  to  elastic  or  elastic-plastic  motions. 

Four  problems  which  involve  the  acceleration  of  metallic  liners 
by  high  explosives  were  solved  numerically  with  the  HEMP  code  and 
these  results  were  compared  with  experimental  observations.  In  all 
cases.  Including  the  case  of  preformed  fragments,  the  results  were  In 
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reasonable  agreement.  In  calculating  the  motion  of  thin  liners,  It  Is 
practical  to  use  only  a  few  mesh  zones  across  the  liner  thickness. 
Therefore,  no  attempt  Is  made  to  resolve  wave  motion  or  stresses  In  the 
metal  liner;  however,  based  upon  the  above  results,  the  gross  motion  of 
the  liner  Is  calculated  reasonably  well. 

From  the  above  comparisons,  we  conclude  that  BRL's  current  ver¬ 
sion  of  the  HEMP  code  Is  operating  correctly  and  that  the  accuracy  Is 
sufficient  to  make  It  useful.  This  code  should  be  applied  to  current 
design  problems  which  Involve  metal  liners  accelerated  by  detonating 
explosives  In  order  to  determine  If  the  number  of  Iterations  required 
by  current  trial  and  error  design  techniques  can  be  significantly  re¬ 
duced. 
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APPENDIX  A 

PROBLEM  FORMULATION  OF  THE  HEMP  CODE 

NOTATION: 

(All  tensors  are  referred  to  a  fixed  Cartesian  coordinate  system.) 

e  specific  Internal  energy 

e. .  components  of  the  rate  of  strain  tensor 

e..  components  of  the  deviator  of  e.. 

*0  tj 

G  shear  modulus 
P  pressure 

S..  components  of  the  deviator  of  0. . 

T'O  IJ 

t  time 

components  of  velocity 

components  of  the  Eulerlan  stress  tensor 

yield  strength  In  simple  tension 

D  material  derivative 
EF 

DS..  deviator  stress  rate  (Jaumann  derivative) 

DT*7 

x  factor  of  proportionality 


p  density 

()p  Implies  plastic  part 
0  Implies  elastic  part 
Kronecker  delta 
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A  description  of  the  problem  formulation  and  numerical  technique 
used  In  the  HEMP  code  can  be  found  in  references  1  and  2.  A  brief  de^ 
scription  is  Included  here  for  easy  reference. 

The  conservation  equations  of  mass,  momentum,  and  energy  may  be 
written  with  reference  to  a  fixed  Cartesian  coordinate  system  as 


fe  +  pu .  .  =  0  j 

Dt 

(A.l) 

(A. 2) 

Oe 

P  Dt  aijeij  , 

(A. 3) 

where  the  standard  convention  of  summation  on  repeated  subscripts  is 
implied  and  a  comma  indicates  partial  differentiation  with  respect  to 
that  space  coordinate.  The  stain  rate  is  defined  as 


13 


*  Kj 


V  ,  .) 


(A. 4) 


For  the  plasticity  theory  used  in  the  code,  the  strain  rate  is  assumed 
to  be  composed  of  an  elastic  and  a  plastic  part  according  to 


The  strain  rate  deviator,  e 


is  defined  by 


e. . 

13 


'13 


3  6ijckk 


(A. 5) 


(A. 6) 


with  similar 'definitions  for  the  elastic,  e.®,  and  plastic,  e?.,  de- 
viators.  The  stress  deviator,  S. .,  and  pressure,  P,  are  given7 by 
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The  von  Mises  yield  function,  f ,  where 


/  = 


(A.8) 


is  used  to  determine  if  the  material  is  in  an  elastic  or  a  plastic  state. 

If  /  <  0  or  /  =  0  and  ^  <  0,  then  the  material  is  in  the  elastic  state  and 
the  above  equations  are  supplemented  by  an  equation  of  state 


P  =  P(p»e)  3 

(A. 9) 

and 

ds  . . 

-M-  =  2Ge® . 

£>t 

O 

H 

•  0> 

* 

(A. 10) 

If  f  =  0  and  ^  >  0,  then  the  material  is  in  the  plastic  state 
preceding  relations  [(A.l)  to  (A.8)]  are  supplemented  by 

and  the 

P  ■  P(p.e)  , 

. 

V  t«7  ’ 

=  °* 

(A. 11) 

The  above  system  of  equations  is  integrated  in  time  by  using  an 
explicit,  finite-difference  technique.  The  derivatives  appearing  in 
the  differential  equations  are  replaced  by  finite-difference  approxi¬ 
mations  as  follows: 


a.  material  derivative 


»  jW^gri^riWWBtHWBWBSSWWW 


where  the  subscript  o  indicates  a  point  fixed  in  the  moving  material, 
b.  partial  space  derivatives 

SI^-TO  +  4<»2-*i)  +  /C(V*2>  ♦  /4<v»,>]  . 

where  the  subscript  notation  is  referred  to  Figure  A- 1<  Quantities 


Figure  A-l.  Finite-Difference  Mesh  Fixed  in  Material 

like  position  and  velocity  are  defined  at  mesh  points,  but  quantities 
like  pressure,  volume,  etc.,  are  defined  at  zone  centers.  In  the  above 
equation,  /  is  a  zone  center  quantity  and  A  is  the  area  of  the  dotted 
polygon.  In  this  formulation,  the  finite-difference  mesh  is  fixed  in 
and  moves  with  the  material.  The  approximation  used  above  is  derived 
from 

¥-sfc*iv- 

A  similar  expression  holds  for  partial  derivatives  with  respect  to  Y. 
c.  stress  rate 


DS  * .  x  i  * 


m&m 
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where  the  prime  indicates  that  ^*is  referred  to  a  Cartesian  coordi¬ 
nate  system  which  translates  and  rotates  with  the  material.  At  time  t, 
the  primed  coordinate  system  is  aligned  with  the  fixed  coordinate  system. 

The  above  finite-difference  equations  are  used  to  reduce  the 
partial  differential  equations  to  algebraic  equations.  Then,  from 
known  conditions  at  one  time  t,  the  variables  can  be  estimated  at 
time  t  +  At.  First,  the  momentum  equations  are  used  to  estimate  the 
velocities  over  the  time  increment  At  from  the  conditions  given  at 
time  t.  The  coordinate  positions  at  t  +  At  for  the  mesh  points  are 
determined  from  these  velocities.  Volumes  at  the  end  of  the  time 
increment  are  determined  from  the  new  coordinate  positions.  The  vis¬ 
cous  pressure  (based  upon  the  von  Neumann-Richtnyer  method)  which  is 
added  to  the  pressure  from  the  equation  of  state  is  calculated  from 


q 


if  *  <  0 


q  =  0  if  iiO, 


or 


where  $  is  the  volumetric  rate  of  change,  A  is  the  zone  area,  and  C  is 
a  constant.  The  strain  rate  is  calculated  from  Equation  (A. 4)  with 
the  known  velocities,  and  the  stress  rate  is  calculated  from  these 
strains  by  Equation  (A. 10).  These  stresses  are  then  used  to  check  the 
yield  condition.  Equation  (A. 8).  If  the  stresses  are  within  the 
elastic  range,  then  these  values  are  correct;  however,  if  the  material 
is  plastic,  these  stresses  are  reduced  by  multiplying  each  component  by 

v^/3  ov(S .  .S . .)  ?  . 

Y'  vj' 

Thus,  the  stress  state  is  reduced  to  fall  on  the  yield  surface,  and 
the  flow  rule  of  Equation  (A.ll)  Is  preserved^  The  density  of 
each  zone  Is  calculated  from  the  continuity  equation  (A.1),  and  e  and 
P  are  obtained  from  the  energy  equation  (A. 3)  and  the  equation  of  state 
(A. 9).  A  time  increment  is  next  calculated  based  upon  the  length  of 
time  required  for  a  dilatational  wave  to  transverse  a  zone.  All 
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quantities  are  now  determined  at  time  t  +  At.  This  procedure  is  con¬ 
tinually  used  to  produce  a  solution  over  the  time  interval  of  interest. 

The  code  can  treat  boundary  conditions  of  either  a  stress  free 
boundary  or  a  specified  velocity  on  the  boundary.  Also,  to  treat 
problems  which  involve  two  materials  which  slide  relative  to  each- 
other,  a  decoupling  of  grid  points  on  the  interface  is  accomplished  in 
a  slide  line  routine. 


